1 Libraries and settings

# ----------------------------------------
# libraries
# ----------------------------------------

library(rtracklayer)
library(GenomicRanges)
library(ggplot2)
library(AnnotationDbi)
library(dplyr)
library(reshape2)
library(UpSetR)
library(GenomicFeatures)
library(kableExtra)
library(knitr)
library(ggrepel)
library(gridExtra)
library(grid)
library(viridis)
library(BindingSiteFinder)
library(ComplexHeatmap)
library(forcats)
library(ggtext)
library(patchwork)
library(tibble)
library(tidyr)
library(dplyr)
library(ggpointdensity)
library(ggsci)
library(ggtext)
library(ggrepel)
library(patchwork)
library(ggrastr)
library(matrixStats)
library(DESeq2)
library(IHW)


source("/Users/melinaklostermann/Documents/projects/AgoCLIP_miR181/mirko_files/mirECLIP/DifferentialBinding/CustomThemes.R")
source("/Users/melinaklostermann/Documents/projects/R_general_functions/theme_paper.R")
source("/Users/melinaklostermann/Documents/projects/AgoCLIP_miR181/mirko_files/mirECLIP/DifferentialBinding/colorPalette.R")

# ----------------------------------------
# settings
# ----------------------------------------

out <- "/Users/melinaklostermann/Documents/projects/AgoCLIP_miR181/R_github/miR181_paper/Figure1/Differential_Binding/"
tpm_cut <- 50
# ----------------------------------------
# ----------------------------------------
# Files
# ----------------------------------------
# ----------------------------------------

# ----------------------------------------
# annotation
# ----------------------------------------

annoDb <- loadDb("/Users/melinaklostermann/Documents/projects/AgoCLIP_miR181/R_github/miR181_paper/Methods/01_Annotation_preprocessing/annotation.db")
gns <- readRDS("/Users/melinaklostermann/Documents/projects/AgoCLIP_miR181/R_github/miR181_paper/Methods/01_Annotation_preprocessing/gene_annotation.rds")

# ----------------------------------------
# AGO binding sites
# ----------------------------------------

bs_wt <- readRDS("/Users/melinaklostermann/Documents/projects/AgoCLIP_miR181/R_github/miR181_paper/Methods/02_AGO_binding_site_definition/2023-04-04-AGO_BS.rds")

bs_ko <- readRDS("/Users/melinaklostermann/Documents/projects/AgoCLIP_miR181/R_github/miR181_paper/Methods/03_KOmir181_AGO_binding_site_definition/2023-04-12-KOmir181_AGO_BS.rds")

# ---------------------------------------- 
# Load clip data
# ----------------------------------------

clipFiles = "/Users/melinaklostermann/Documents/projects/AgoCLIP_miR181/pipe_output_22_02_14/non-chimeric/bw"
clipFiles = list.files(clipFiles, pattern = ".bw$", full.names = TRUE)
clipFiles = clipFiles[!grepl("Inp", clipFiles)]
clipFiles = clipFiles[!grepl("miR181_", clipFiles)]
clipFilesP = clipFiles[grep(clipFiles, pattern = "plus")]
clipFilesM = clipFiles[grep(clipFiles, pattern = "minus")]

2 What was done?

3 Combine binding sites

Binding sites from both conditions do either overlap exactly, overlap partially or don’t overlap at all. In the following partial overlaps are resolved by re-centering the binding sites based on the highest crosslink signal of both conditions.

# merge BS from both conditions
bsMerge = unique(c(resize(bs_wt, fix = "center", width = 1),
                   resize(bs_ko, fix = "center", width = 1)))

# Organize clip data in dataframe
colData = data.frame(
  id = c(1:6),
  condition = factor(c(rep("KO",3), rep("WT",3)),
                     levels = c("KO", "WT")),
                     clPlus = clipFilesP,
                     clMinus = clipFilesM)

# Make BindingSiteFinder object
bds = BSFDataSetFromBigWig(ranges = bsMerge, meta = colData)


# Make unified binding sites from both conditions
bdsMerge <- makeBindingSites(object = bds, bsSize = 7, minWidth = 0,
                        minCrosslinks = 0, minClSites = 0, centerIsSummit = FALSE)

df = getSummary(bdsMerge)
kable(myNumberFormat(df), caption = "Merge and combine") 
Table 1: Merge and combine
Option nRanges
inputRanges 37,575
mergeCrosslinkSites 28,346
minCrosslinks 28,346
minClSites 28,346
centerIsClSite 28,253
centerIsSummit NA

Combination of both sets of binding sites results in a single set (N=28,253), where each binding site was either seen in WT, KO or both conditions. The figure below shows the set sizes in more detail.

Annotate combined binding sites with PureCLIP score, WT and KO status information.

# annotate with pureclip score
rng = getRanges(bds)
rng$additionalScores = NULL
mcols(rng)$score = rng$scoreSum
bdsFinal = annotateWithScore(bdsMerge, rng)

# annotate with condition support
rngFinal = getRanges(bdsFinal)
r1 = resize(bs_wt, fix = "center", width = 1)
r2 = resize(bs_ko, fix = "center", width = 1)
condition_support = data.frame(WT = countOverlaps(rngFinal, r1),
                KO = countOverlaps(rngFinal, r2))
mcols(rngFinal) = cbind(mcols(rngFinal), condition_support)

# transfer the geneID 
olsWT = findOverlaps(rngFinal, bs_wt)
olsKO = findOverlaps(rngFinal, bs_ko)
rngFinal$geneID = NA
rngFinal$geneID[queryHits(olsWT)] = bs_wt$geneID[subjectHits(olsWT)]
rngFinal$geneID[queryHits(olsKO)] = bs_ko$geneID[subjectHits(olsKO)]

# transfer gene names
rngFinal$geneName = NA
rngFinal$geneName[queryHits(olsWT)] = bs_wt$geneName[subjectHits(olsWT)]
rngFinal$geneName[queryHits(olsKO)] = bs_ko$geneName[subjectHits(olsKO)]

# transfer transcript region
rngFinal$region = NA
rngFinal$region[queryHits(olsWT)] = bs_wt$region[subjectHits(olsWT)]
rngFinal$region[queryHits(olsKO)] = bs_ko$region[subjectHits(olsKO)]

# set final range 
bdsFinal = setRanges(bdsFinal, rngFinal)

4 Perform Differntial binding

In order to perform differential binding analysis, we use DEseq2. The design formula contains both the number of crosslinks per binding site and the background crosslinks in the whole gene of the binding site.

4.2 Differential analysis of binding sites

### ============================================================================
### Run test with DESeq
### ----------------------------------------------------------------------------
#
count_matrix = as.data.frame(mcols(countObj))
count_matrix = count_matrix %>% select(starts_with("counts"))

# internal modification for col data
colDataMod = rbind.data.frame(colData,colData)
colDataMod$type = c(rep("bs", nrow(colData)),
                    rep("bg", nrow(colData)))
colDataMod$condition = factor(colDataMod$condition, levels = c("KO", "WT"))
colDataMod$type = factor(colDataMod$type, levels = c("bs", "bg"))

# create SE object
se = SummarizedExperiment(assays = list(counts = as.matrix(count_matrix)),
                          rowRanges = granges(countObj), colData = colDataMod)

# set design - YOUs design
ddsBs = DESeqDataSet(se, design = ~condition + type + condition:type)
ddsBs$condition = relevel(ddsBs$condition, "WT")
ddsBs$type = relevel(ddsBs$type, "bg")
ddsBs = DESeq(ddsBs, test="LRT", reduced = ~ condition + type)
resBs = results(ddsBs, name = "conditionKO.typebs") # needs relevel of type to bg
resShrinkBs = lfcShrink(ddsBs, res = resBs, coef = "conditionKO.typebs", type = "ashr")


# deseq fetch results 
diff_bs = countObj
colnames(resBs) = paste0("resBs.", colnames(resBs))
idx = match(names(diff_bs), rownames(resBs))
mcols(diff_bs) = cbind(mcols(diff_bs), resBs[idx,])

4.3 Differential analysis of background

### ============================================================================
### Run Background test with DESeq
### ----------------------------------------------------------------------------
#
# construct background signal dataframe to test for the background level change
count_matrix_bg = as.data.frame(mcols(countObj)) %>% 
    select(starts_with("counts.bg")) %>%
    unique() 
#rownames(count_matrix_bg) = sapply(strsplit(rownames(count_matrix_bg),"\\."), `[`, 1)

# create SE object
seBg = SummarizedExperiment(assays = list(counts = as.matrix(count_matrix_bg)), colData = colData)

# DESeq design
ddsBg = DESeqDataSet(seBg, design = ~ condition)
ddsBg$condition = relevel(ddsBg$condition, "WT")
ddsBg = DESeq(ddsBg)
resBg = results(ddsBg, contrast = c("condition", "KO", "WT"))
resShrinkBg = lfcShrink(ddsBg, res = resBg, contrast = c("condition", "KO", "WT"), type = "ashr")

# add results to final differential object
colnames(resBg) = paste0("resBg.", colnames(resBg))
idx = match(names(diff_bs), rownames(resBg))
mcols(diff_bs) = cbind(mcols(diff_bs), resBg[idx,])
# ### ============================================================================
# ### Run enhanced background test with DESeq
# ### ----------------------------------------------------------------------------
# #
# # load WT
# peaksInitialWT = "/Users/mirko/Projects/mirEClip/01_eCLIP/pureclip/IP_WT_pureclip_sites.bed"
# peaksInitialWT = import(con = peaksInitialWT, format = "BED", extraCols=c("additionalScores" = "character"))
# peaksInitialWT = keepStandardChromosomes(peaksInitialWT, pruning.mode = "coarse")
# # load KO
# peaksInitialKO = "/Users/mirko/Projects/mirEClip/01_eCLIP/pureclip/IP_KO_pureclip_sites.bed"
# peaksInitialKO = import(con = peaksInitialKO, format = "BED", extraCols=c("additionalScores" = "character"))
# peaksInitialKO = keepStandardChromosomes(peaksInitialKO, pruning.mode = "coarse")
# # peaks all
# peaksAll = unique(c(peaksInitialKO, peaksInitialWT)) 
# bdsPeaksAll = setRanges(bds, peaksAll)
# 
# countObj2 = makeBsBackgroundMatrix(object = bdsPeaksAll, geneRanges = gns, offset = 5, matchBy = "range") 
# # construct background signal dataframe to test for the background level change
# mAllPeaks = as.data.frame(mcols(countObj2))
# mAllPeaks = mAllPeaks %>% 
#     select(starts_with("counts.bg")) %>%
#     unique()
# rownames(mAllPeaks) = sapply(strsplit(rownames(mAllPeaks),"\\."), `[`, 1)
# 
# # create SE object
# seBg2 = SummarizedExperiment(assays = list(counts = as.matrix(mAllPeaks)), colData = colData)
# 
# # DESeq design
# ddsBg2 = DESeqDataSet(seBg2, design = ~ condition)
# ddsBg2$condition = relevel(ddsBg2$condition, "WT")
# ddsBg2 = DESeq(ddsBg2)
# resBg2 = results(ddsBg2, contrast = c("condition", "KO", "WT"))
# resBg2 = lfcShrink(ddsBg2, res = resBg2, contrast = c("condition", "KO", "WT"), type = "ashr")


### ============================================================================
### Run Gene test with DESeq
### ----------------------------------------------------------------------------
#
# # construct gene signal dataframe to test for total gene level change
# m = mcols(countObj) %>%
#     as.data.frame() %>%
#     select(starts_with("counts.bs"), geneID) %>%
#     group_by(geneID) %>%
#     summarise_all(sum) %>%
#     tibble::column_to_rownames("geneID")
# 
# # create SE object
# seGene = SummarizedExperiment(assays = list(counts = as.matrix(m)), colData = colData)
# 
# # DESeq design
# ddsGene = DESeqDataSet(seGene, design = ~ condition)
# ddsGene$condition = relevel(ddsGene$condition, "WT")
# ddsGene = DESeq(ddsGene)
# resGene = results(ddsGene, contrast = c("condition", "KO", "WT"))
# resGene = lfcShrink(ddsGene, res = resGene, contrast = c("condition", "KO", "WT"), type = "ashr")


### ============================================================================
### Export results
### ----------------------------------------------------------------------------
#

4.4 Initial results

### ============================================================================
### Numbers overview table
### ----------------------------------------------------------------------------
### 
df = data.frame(Name = c("Tested", "Not Significant", "Significant", "Sig+Up", "Sig+Down"),
                N = c(nrow(resBs),
                      nrow(subset(resBs, resBs.padj >= 0.05)),
                      nrow(subset(resBs, resBs.padj < 0.05)),
                      nrow(subset(resBs, resBs.padj < 0.05 & resBs.log2FoldChange > 0)),
                      nrow(subset(resBs, resBs.padj < 0.05 & resBs.log2FoldChange < 0))
                      ))
df$Per = c(100,
           nrow(subset(resBs, resBs.padj >= 0.05)) / nrow(resBs) * 100,
           nrow(subset(resBs, resBs.padj < 0.05)) / nrow(resBs) * 100,
           nrow(subset(resBs, resBs.padj < 0.05 & resBs.log2FoldChange > 0)) / nrow(subset(resBs, resBs.padj < 0.05)) * 100,
           nrow(subset(resBs, resBs.padj < 0.05 & resBs.log2FoldChange < 0)) / nrow(subset(resBs, resBs.padj < 0.05)) * 100
           )
df$Per = round(df$Per, digits = 2)
colnames(df) = c("Name", "Count (#N)", "Percentage (%)")

kable(myNumberFormat(df), caption = "Results by numbers. Significant based on IHW adjusted P value threshold 0.05.") 
Table 2: Results by numbers. Significant based on IHW adjusted P value threshold 0.05.
Name Count (#N) Percentage (%)
Tested 28,253 100.00
Not Significant 27,115 95.97
Significant 1,138 4.03
Sig+Up 296 26.01
Sig+Down 842 73.99

4.5 Gene expression filter

Some genes might not be expressed in both conditions. Binding sites on a gene that is for example only expressed in the WT but not in the KO condition will all appear to be down regulated. In reality we can not tell if such sites changed, because the hosting gene is not expressed. To prevent these effects I filtered all genes with a merged binding site (so the combined set from both conditions) before running the differential binding search. As cutoff a TPM > 50 is used.

# cutoff based on gene

counts = counts(ddsBg, normalized = FALSE)
# Extract all exons of a gen
exonsByGene <- exonsBy(annoDb, by = "gene")
# reduce overlapping exons to a single region
reducedExonsByGene <- reduce(exonsByGene)
# Approximate the gene length as the sum of the its exons lengths
geneLengths <- sum(width(reducedExonsByGene))
# Adapt gene IDs
names(geneLengths) = sapply(strsplit(names(geneLengths),"\\."), `[`, 1)
# Re-order the vector of gene lengths to match the order in the counts
counts_genes <- sapply(strsplit(rownames(counts),"\\."),  `[`, 1)
geneLengths <- geneLengths[match(counts_genes, names(geneLengths))]
# First step1 in TPM calculation
tpms <- counts / geneLengths
# Second step in TPM calculation
tpms <- t(t(tpms) * 1e6 / colSums(tpms, na.rm = T)) %>% as.data.frame()
# # Adapt gene IDs
rownames(tpms) = sapply(strsplit(rownames(tpms),"\\."), `[`, 1)
colnames(tpms) = paste0("tpm.", colnames(tpms) )

# add tpm to final object
idx <- match(diff_bs$geneID, rownames(tpms))
mcols(diff_bs) <- cbind(mcols(diff_bs), tpms[idx,])

df = tpms %>%
    as.data.frame() %>%
    tibble::rownames_to_column("geneID") %>%
    pivot_longer(-geneID) %>%
    group_by(name) %>% 
    mutate(rank = rank(value, ties.method = "first"))


ggplot(df, aes(x = rank, y = value)) +
    ggrastr::rasterise(geom_pointdensity(size = 2), dpi = 300) +
    facet_wrap(~name) +
    theme_nice() +
    geom_hline(yintercept = 50, linetype = "dashed") +
    scale_y_log10() +
    scale_color_viridis(option = "rocket") +
    theme(legend.position = "top") +
    guides(color = guide_colorbar(title.position = 'top', title.hjust = 0.5,
                                  barwidth = unit(20, 'lines'), barheight = unit(.5, 'lines'))) +
    labs(
        title = "TPMs per gene per sample",
        x = "Rank",
        y = "TPMs",
        color = "#Points"
        )
TPMs per gene for each sample, sorted by rank on log10 scale. Cutoff is TPM > 50

Figure 1: TPMs per gene for each sample, sorted by rank on log10 scale
Cutoff is TPM > 50

A gene is considered expressed by a condition if at least two of the three replicates have a TPM over 50. Each gene must be found expressed in both conditions to be considered for the differential binding analysis.

# filter for tpm cutoff in 2 samples per condition
diff_bs$BS_ID <- names(diff_bs)
diff_bs <- as.data.frame(diff_bs) %>%
  rowwise() %>%
  mutate(tpm_support_KO = sum(c((tpm.counts.bg.1_KO > tpm_cut),
                             (tpm.counts.bg.2_KO > tpm_cut),
                             (tpm.counts.bg.3_KO > tpm_cut))),
         tpm_support_WT = sum(c((tpm.counts.bg.4_WT > tpm_cut),
                             (tpm.counts.bg.5_WT > tpm_cut),
                             (tpm.counts.bg.6_WT > tpm_cut))),
         tpm_supported = (tpm_support_KO > 1) & (tpm_support_WT > 1)) 
                           

# make upset plot
genesExpInWT = diff_bs %>% subset(tpm_support_WT > 1) %>% pull(geneID) %>% unique()
genesExpInKO = diff_bs %>% subset(tpm_support_KO > 1) %>% pull(geneID) %>% unique()
  
l = list(genesExpInWT = genesExpInWT, genesExpInKO = genesExpInKO)
m = make_comb_mat(l)

ha = HeatmapAnnotation(
  "Intersections" = anno_barplot(comb_size(m), border = FALSE, gp = gpar(fill = "#a6a6a6"), height = unit(6, "cm"))
)
ha2 = HeatmapAnnotation(
  "set sizes" = anno_barplot(set_size(m), border = FALSE, gp = gpar(fill = "#a6a6a6"), width = unit(4, "cm")),
  which = "row"
)

ht = UpSet(m,
           comb_order = order(comb_size(m), decreasing = T),
           top_annotation = ha,
           right_annotation = ha2,
           comb_col = "#003399", bg_col = "white", pt_size = unit(.5, "cm") ,
           border = T, lwd = 2, bg_pt_col = "#a6a6a6"
)

ss = set_size(m)
cs = comb_size(m)
ht = draw(ht,  padding = unit(c(0, 0, 10, 0), "mm"))
od = column_order(ht)
decorate_annotation("Intersections", {
    grid.text(format(cs[od], big.mark = ".", decimal.mark = ","), x = seq_along(cs), y = unit(cs[od], "native") + unit(.1, "pt"), 
        default.units = "native", just = c("left", "bottom"), 
        gp = gpar(fontsize = 8, col = "black"), rot = 45)
})
Intersection of genes expressed in both conditions.

Figure 2: Intersection of genes expressed in both conditions

# subset by cutoff
diff_bs <- diff_bs %>% subset(tpm_supported )

The expression filtering resulted in 3,779 target genes, with 26,462 binding sites.

5 Results (after Gene Expression filter)

### ============================================================================
### Numbers overview table
### ----------------------------------------------------------------------------
### 
df = data.frame(Name = c("Tested", "Not Significant", "Significant", "Sig+Up", "Sig+Down"),
                N = c(nrow(diff_bs),
                      nrow(subset(diff_bs, resBs.padj >= 0.05)),
                      nrow(subset(diff_bs, resBs.padj < 0.05)),
                      nrow(subset(diff_bs, resBs.padj < 0.05 & resBs.log2FoldChange > 0)),
                      nrow(subset(diff_bs, resBs.padj < 0.05 & resBs.log2FoldChange < 0))
                      ))
df$Per = c(100,
           nrow(subset(diff_bs, resBs.padj >= 0.05)) / nrow(diff_bs) * 100,
           nrow(subset(diff_bs, resBs.padj < 0.05)) / nrow(diff_bs) * 100,
           nrow(subset(diff_bs, resBs.padj < 0.05 & resBs.log2FoldChange > 0)) / nrow(subset(diff_bs, resBs.padj < 0.05)) * 100,
           nrow(subset(diff_bs, resBs.padj < 0.05 & resBs.log2FoldChange < 0)) / nrow(subset(diff_bs, resBs.padj < 0.05)) * 100
           )
df$Per = round(df$Per, digits = 2)
colnames(df) = c("Name", "Count (#N)", "Percentage (%)")

kable(myNumberFormat(df), caption = "Results by numbers. Significant based on IHW adjusted P value threshold 0.05.") 
Table 3: Results by numbers. Significant based on IHW adjusted P value threshold 0.05.
Name Count (#N) Percentage (%)
Tested 26,462 100.00
Not Significant 25,388 95.94
Significant 1,074 4.06
Sig+Up 278 25.88
Sig+Down 796 74.12

5.1 Decide for P value cutoff

d = lapply(seq(0, 1, by = 0.05), function(cutoff){
    diff_bs %>%
        select(region, resBs.log2FoldChange, resBs.padj) %>%
        rename("lfc" = "resBs.log2FoldChange", "padj" = "resBs.padj") %>%
        subset(padj <= cutoff) %>%
        mutate(dir = ifelse(lfc > 0, "Up", "Down")) %>%
        select(region, dir) %>%
        mutate(cutoff = cutoff)
})
d = do.call("rbind", d)

df = d %>%
    group_by(region, dir, cutoff) %>% 
    tally() %>%
    ungroup() %>%
    group_by(cutoff, region) %>%
    mutate(sum = sum(n)) %>%
    arrange(cutoff) %>%
    mutate(percentage = round(n/sum, digits = 4)*100) %>%
    mutate(dir = factor(dir, level = c("Up", "Down"))) %>%
    arrange(dir)

p1 = ggplot(df, aes(x = cutoff, y = percentage, color = dir, fill = dir, group = dir)) +
    geom_line() +
    geom_point(size = 4, shape = 21) +
    theme_pub() +
    scale_fill_nice_ud_b() +
    scale_color_nice_ud_d() +
    theme(legend.position = "top") +
    labs(
        title = "Regulation enrichment",
        x = "Adjusted P value cutoffs",
        y = "Percentage of significantly regulated binding sites",
        color = "Direction",
        fill = "Direction"
    ) +
    theme(axis.text.x=element_text(angle = 45, hjust = 1, vjust = 1)) +
    facet_wrap(~region) +
    geom_vline(xintercept = 0.075, linetype = "dashed")

p2 = ggplot(df, aes(x = cutoff, y = n, fill = dir, color = dir, group = dir)) +
    geom_col() +
    theme_pub() +
    scale_fill_nice_ud_b() +
    scale_color_nice_ud_d() +
    theme(legend.position = "top") +
    theme(axis.text.x=element_text(angle = 45, hjust = 1, vjust = 1)) +
    facet_wrap(~region, scales = "free_y") +
    geom_vline(xintercept = 0.075, linetype = "dashed") +
    labs(
        title = "Number of binding site per adjusted P value cutoff",
        x = "Adjusted P value cutoffs",
        y = "Count",
        color = "Direction",
        fill = "Direction"
    ) 

p1 / p2
Regulation enrichment per adjusted P value cutoff. The dashed line indicates a cutoff of 0.05.

Figure 3: Regulation enrichment per adjusted P value cutoff
The dashed line indicates a cutoff of 0.05.

# P value hist
p1 = ggplot(diff_bs, aes(x = resBs.pvalue)) +
    geom_histogram(bins = 50, fill = "#999999", color = "#4d4d4d") +
    theme_nice() +
    labs(
        x = "P value",
        y = "Count"
    )

p2 = ggplot(diff_bs, aes(x = resBs.padj)) +
    geom_histogram(bins = 50, fill = "#999999", color = "#4d4d4d") +
    theme_nice() +
    geom_vline(xintercept = 0.05) +
    labs(
        x = "Adj P value",
        y = "Count"
    )

p1 + p2
P value distributions.

Figure 4: P value distributions

df2 = d %>%
    subset(., cutoff == 0.05) %>%
    subset(., !is.na(region)) %>%
    subset(., region != "outside")%>%
    group_by(region, dir) %>% 
    tally() %>%
    ungroup() %>%
    mutate(sum = sum(n)) %>%
    mutate(percentage = round(n/sum, digits = 4)*100) 
df2 = df2 %>%
    mutate(region = factor(region, levels = c("utr5", "intron", "cds", "utr3"))) %>%
    arrange(region) 

p1 = ggplot(df2, aes(x = dir, y = percentage, fill = region)) + 
    geom_col(position = "stack") +
    theme_pub() +
    scale_fill_npg() +
    labs(
        title = "Significantly regulated binding sites per region and direction",
        x = "Direction of regulation",
        y = "Percentage",
        fill = "Region"
    ) + 
    coord_flip() +
    theme(legend.position = "top") 

p2 = ggplot(df2, aes(x = dir, y = percentage, fill = region)) + 
    geom_col(position = "fill", lwd = 2) +
    theme_pub() +
    scale_fill_npg() +
    labs(
        title = "Significantly regulated binding sites per region and direction",
        x = "Direction of regulation",
        y = "Percentage",
        fill = "Region"
    ) + 
    coord_flip() +
    theme(legend.position = "top") 
    # annotate(geom = "rect", xmin = 0.55, xmax = 1.45, ymin = 0, ymax = 1, fill = "NA", color = "#874C62", lwd = 2) +
    # annotate(geom = "rect", xmin = 1.55, xmax = 2.45, ymin = 0, ymax = 1, fill = "NA", color = "#68b1a5", lwd = 2)

p1 + p2
Significantly regulated binding sites per region and direction. The percentage is calculated from all significant binding sites at adusted P value cutoff of 0.05.

Figure 5: Significantly regulated binding sites per region and direction
The percentage is calculated from all significant binding sites at adusted P value cutoff of 0.05.

5.1.1 Decide for fold-change cutoff

d = lapply(log2(seq(1, 10, by = 0.5)), function(cutoff){
    diff_bs %>%
        filter(resBs.padj < 0.05 & abs(resBs.log2FoldChange) > cutoff) %>%
        mutate(dir = factor(ifelse(resBs.log2FoldChange > 0, "Up", "Down"), level = c("Up", "Down"))) %>%
        select(dir) %>%
        group_by(dir) %>%
        tally() %>%
        # mutate(cutoff = paste0("log2(", 2^cutoff, ")")) %>%
        mutate(cutoff = 2^cutoff ) %>%
        # mutate(cutoff = factor(cutoff, levels = c(paste0("log2(", seq(1, 10, by = 0.5), ")")))) %>%
        mutate(cutoff = factor(cutoff, levels = c(paste0(seq(1, 10, by = 0.5))))) %>%
        arrange(cutoff)
})
d = do.call("rbind", d)
df = d %>%
    group_by(cutoff) %>%
    mutate(sum = sum(n)) %>%
    arrange(cutoff) %>%
    mutate(percentage = round(n/sum, digits = 4)*100) 

p1 = ggplot(df, aes(x = cutoff, y = percentage, color = dir, fill = dir, group = dir)) +
    geom_line() +
    geom_point(size = 4, shape = 21) +
    theme_pub() +
    scale_fill_nice_ud_b() +
    scale_color_nice_ud_d() +
    theme(legend.position = "top") +
    # theme(axis.text.x=element_text(angle = 45, hjust = 1, vjust = 1)) +
    labs(
        title = "Regulation enrichment",
        x = "Fold-change cutoffs (log2)",
        y = "Percentage of regulated binding sites",
        color = "Direction",
        fill = "Direction"
    ) 

p2 = ggplot(df, aes(x = cutoff, y = n, fill = dir, color = dir, group = dir)) +
    geom_col() +
    theme_pub() +
    scale_fill_nice_ud_b() +
    scale_color_nice_ud_d() +
    theme(legend.position = "top") +
    # theme(axis.text.x=element_text(angle = 45, hjust = 1, vjust = 1)) +
    labs(
        title = "Number of binding site per fold-change cutoff",
        x = "Fold-change cutoffs (log2)",
        y = "Count",
        color = "Direction",
        fill = "Direction"
    ) 
    
p1 + p2 
Enrichment of the direction of regulation. Based on differnt cutoffs for the fold-change per binding site.

Figure 6: Enrichment of the direction of regulation
Based on differnt cutoffs for the fold-change per binding site.

df = diff_bs %>%
    mutate(dir = factor(ifelse(resBs.log2FoldChange > 0, "Up", "Down"), level = c("Up", "Down"))) %>%
    arrange(dir)
p1 = ggplot(df, aes(x = resBs.log2FoldChange)) +
    geom_histogram(bins = 50, fill = "#999999", color = "#4d4d4d") +
    theme_nice() +
    labs(
        x = "Fold-change (log2)",
        y = "Count"
    )

df = diff_bs %>%
    mutate(sig = factor(ifelse(resBs.padj < 0.05, TRUE, FALSE))) %>%
    mutate(dir = factor(ifelse(resBs.log2FoldChange > 0, "Up", "Down"), level = c("Up", "Down"))) %>%
    arrange(dir)
p2 = ggplot(subset(df, sig == TRUE), aes(x = resBs.log2FoldChange, fill = dir, color = dir)) +
    geom_histogram(bins = 100) +
    scale_fill_nice_ud_b() +
    scale_color_nice_ud_d() +
    theme_nice() +
    theme(legend.position = "top") +
    labs(
        x = "Fold-change (log2)",
        y = "Count",
        color = "Regulation",
        fill = "Regulation"
    ) 

df = diff_bs %>%
    mutate(sig = factor(ifelse(resBs.padj < 0.05 & abs(resBs.log2FoldChange) > log2(2), TRUE, FALSE))) %>%
    mutate(dir = factor(ifelse(resBs.log2FoldChange > 0, "Up", "Down"), level = c("Up", "Down"))) %>%
    arrange(dir)
p3 = ggplot(subset(df, sig == TRUE), aes(x = resBs.log2FoldChange, fill = dir, color = dir)) +
    geom_histogram(bins = 100) +
    scale_fill_nice_ud_b() +
    scale_color_nice_ud_d() +
    theme_nice() +
    theme(legend.position = "top") +
    labs(
        x = "Fold-change (log2)",
        y = "Count",
        color = "Regulation",
        fill = "Regulation"
    ) 
p1 + p2 + p3
Fold-change distributions. Fold-Change cutoff in the last plot = 2).

Figure 7: Fold-change distributions
Fold-Change cutoff in the last plot = 2).

5.1.2 Final results plots

# MA + Volcano plots
df = diff_bs %>%
    mutate(sig = (ifelse(resBs.padj < 0.05 & abs(resBs.log2FoldChange) > log2(2), TRUE, FALSE))) %>%
    mutate(dir = factor(ifelse(resBs.log2FoldChange > 0, "Up", "Down"), level = c("Up", "Down"))) %>%
    mutate(sigDir = ifelse(sig == TRUE & dir == "Up", "Up", ifelse(sig == TRUE & dir == "Down", "Down", "Not"))) %>%
    mutate(sigDir = factor(sigDir, levels = c("Not", "Up", "Down"))) %>%
    arrange(sigDir)
p1 = ggplot(df, aes(x = log2(resBs.baseMean), y = resBs.log2FoldChange, color = sigDir, fill = sigDir)) +
  geom_point_rast(shape = 21, stroke = 0.5, size = 1.5) +
  scale_fill_nice_udn_b() +
  scale_color_nice_udn_d() +
  geom_hline(yintercept = 0, color = "black", alpha = .5) +
  theme_nice() +
  guides(color = guide_legend(override.aes = list(size = 4))) + 
  theme(legend.key.size = unit(1, 'cm'), legend.position = "top") +
  labs(
      x = "Base mean", 
      y = "Fold-change (log2)", 
      color = "Regulation",
      fill = "Regulation")

p2 = ggplot(df, aes(x = resBs.log2FoldChange, y = -log10(resBs.padj), color = sigDir, fill = sigDir)) +
    geom_point_rast(shape = 21, stroke = 0.5, size = 1.5) +
    scale_fill_nice_udn_b() +
    scale_color_nice_udn_d() +
    geom_vline(xintercept = 0, color = "black", alpha = .5) +
    theme_nice() +
    guides(color = guide_legend(override.aes = list(size = 4))) + 
    theme(legend.key.size = unit(1, 'cm'), legend.position = "top") +
    labs(
        x = "Fold-change (log2)", 
        y = "Adjusted P value (-log10)", 
        color = "Regulation",
        fill = "Regulation")

(p1 + p2) + plot_layout(guides = "collect") & theme(legend.position = "top")
MA and Volcano plots. With adjusted P value and fold-change cutoffs. Simple version.

Figure 8: MA and Volcano plots
With adjusted P value and fold-change cutoffs. Simple version.

p2 <- p2+theme_paper()
ggsave(p2, file= paste0(out, "Figure1X_Differntial_binding_vulcano.pdf"), height = 6, width = 6, units = "cm")
# MA + Volcano plots
df = diff_bs %>%
    mutate(sig = (ifelse(resBs.padj < 0.05 & abs(resBs.log2FoldChange) > log2(2), TRUE, FALSE))) %>%
    mutate(dir = factor(ifelse(resBs.log2FoldChange > 0, "Up", "Down"), level = c("Up", "Down"))) %>%
    mutate(sigDir = ifelse(sig == TRUE & dir == "Up", "Up", ifelse(sig == TRUE & dir == "Down", "Down", "Not"))) %>%
    mutate(sigDir = factor(sigDir, levels = c("Not", "Up", "Down"))) %>%
    arrange(sigDir)

dfN = df %>% subset(sig == TRUE) %>% select(dir) %>% group_by(dir) %>% tally()

dfLabMa_up = df %>% subset(sigDir == "Up") %>% 
    arrange(desc(resBs.log2FoldChange) ) %>% 
    head(5)
dfLabMa_down = df %>% subset(sigDir == "Down") %>% 
    arrange((resBs.log2FoldChange) ) %>% 
    head(5)


p1 = ggplot(df, aes(x = log2(resBs.baseMean), y = resBs.log2FoldChange, color = sigDir, fill = sigDir)) +
    geom_point_rast(shape = 21, stroke = 0.5, size = 1.5) +
    scale_fill_nice_udn_b() +
    scale_color_nice_udn_d() +
    geom_hline(yintercept = 0, color = "black", alpha = .5) +
    theme_nice() +
    guides(color = guide_legend(override.aes = list(size = 4))) + 
    theme(legend.key.size = unit(1, 'cm'), legend.position = "top") +
    labs(
        x = "Base mean", 
        y = "Fold-change (log2)", 
        color = "Regulation",
        fill = "Regulation") +
    geom_text_repel(data = dfLabMa_up, aes(x = log2(resBs.baseMean), y = resBs.log2FoldChange, color = sigDir, fill = sigDir, label = paste0(geneName,"-", sapply(strsplit(BS_ID,"\\."), `[`, 2))),
                    min.segment.length = 0, nudge_y = 2, size = 2.5) +
    geom_text_repel(data = dfLabMa_down, aes(x = log2(resBs.baseMean), y = resBs.log2FoldChange, color = sigDir, fill = sigDir, label = paste0(geneName,"-", sapply(strsplit(BS_ID,"\\."), `[`, 2))),
                    min.segment.length = 0, nudge_y = -2, size = 2.5) +
    ylim(-8,8) +
    annotate(geom = "text", label = paste0("#N=",dfN[1,][[2]]), x = 12.5, y = 5) +
    annotate(geom = "text", label = paste0("#N=",dfN[2,][[2]]), x = 12.5, y = -5)

dfLabVol_up = df %>% filter(sigDir == "Up") %>% 
    arrange((resBs.padj) ) %>% 
    head(5)
dfLabVol_down = df %>% filter(sigDir == "Down") %>% 
    arrange((resBs.padj) ) %>% 
    head(5)  

p2 = ggplot(df, aes(x = resBs.log2FoldChange, y = -log10(resBs.padj), color = sigDir, fill = sigDir)) +
    geom_point_rast(shape = 21, stroke = 0.5, size = 1.5) +
    scale_fill_nice_udn_b() +
    scale_color_nice_udn_d() +
    geom_vline(xintercept = 0, color = "black", alpha = .5) +
    theme_nice() +
    guides(color = guide_legend(override.aes = list(size = 4))) + 
    theme(legend.key.size = unit(1, 'cm'), legend.position = "top") +
    labs(
        x = "Fold-change (log2)", 
        y = "Adjusted P value (-log10)", 
        color = "Regulation",
        fill = "Regulation") +
    geom_text_repel(data = dfLabVol_up, aes(x = resBs.log2FoldChange, y = -log10(resBs.padj), color = sigDir, fill = sigDir, label = paste0(geneName,"-", sapply(strsplit(BS_ID,"\\."), `[`, 2))),
                    min.segment.length = 0, nudge_x = 2, size = 2.5) +
    geom_text_repel(data = dfLabVol_down, aes(x = resBs.log2FoldChange, y = -log10(resBs.padj), color = sigDir, fill = sigDir, label = paste0(geneName,"-", sapply(strsplit(BS_ID,"\\."), `[`, 2))),
                    min.segment.length = 0, nudge_x = -2, size = 2.5) +
    xlim(-6,6) +
    annotate(geom = "text", label = paste0("#N=",dfN[1,][[2]]), y = 50, x = 2.5) +
    annotate(geom = "text", label = paste0("#N=",dfN[2,][[2]]), y = 50, x = -2.5)

(p1 + p2) + plot_layout(guides = "collect") & theme(legend.position = "top")
MA and Volcano plots. With adjusted P value and fold-change cutoffs. Additional annotaitons.

Figure 9: MA and Volcano plots
With adjusted P value and fold-change cutoffs. Additional annotaitons.

# MA + Volcano plots
df = diff_bs %>%
    mutate(sig = (ifelse(resBs.padj < 0.05 & abs(resBs.log2FoldChange) > log2(2), TRUE, FALSE))) %>%
    mutate(dir = factor(ifelse(resBs.log2FoldChange > 0, "Up", "Down"), level = c("Up", "Down"))) %>%
    mutate(sigDir = ifelse(sig == TRUE & dir == "Up", "Up", ifelse(sig == TRUE & dir == "Down", "Down", "Not"))) %>%
    mutate(sigDir = factor(sigDir, levels = c("Not", "Up", "Down"))) %>%
    arrange(sigDir)

dfN = df %>% filter(sig == TRUE) %>% select(dir) %>% group_by(dir) %>% tally()

selectedGenes = c("Zfp36l1","Zfp36l2","Dusp10","Ddit4","Gnb4","S1pr1")

dfLabMa_up = df %>% filter(geneName %in% selectedGenes & resBs.log2FoldChange > 0 & sig == TRUE)
dfLabMa_down = df %>% filter(geneName %in% selectedGenes & resBs.log2FoldChange < 0 & sig == TRUE)

p1 = ggplot(df, aes(x = log2(resBs.baseMean), y = resBs.log2FoldChange, color = sigDir, fill = sigDir)) +
    geom_point_rast(shape = 21, stroke = 0.5, size = 1.5) +
    scale_fill_nice_udn_b() +
    scale_color_nice_udn_d() +
    geom_hline(yintercept = 0, color = "black", alpha = .5) +
    theme_nice() +
    guides(color = guide_legend(override.aes = list(size = 4))) + 
    theme(legend.key.size = unit(1, 'cm'), legend.position = "top") +
    labs(
        x = "Base mean", 
        y = "Fold-change (log2)", 
        color = "Regulation",
        fill = "Regulation") +
    # geom_text_repel(data = dfLabMa_up, aes(x = log2(res.baseMean), y = res.log2FoldChange, color = sigDir, fill = sigDir, label = paste0(geneName,"-", sapply(strsplit(BsID,"\\."), `[`, 2))),
    #                 min.segment.length = 0, nudge_y = 6, size = 2.5) +
    geom_text_repel(data = dfLabMa_down, aes(x = log2(resBs.baseMean), y = resBs.log2FoldChange, color = sigDir, fill = sigDir, label = paste0(geneName,"-", sapply(strsplit(BS_ID,"\\."), `[`, 2))),
                    min.segment.length = 0, nudge_y = -6, size = 2.5) +
    ylim(-7,7) +
    annotate(geom = "text", label = paste0("#N=",dfN[1,][[2]]), x = 12.5, y = 5) +
    annotate(geom = "text", label = paste0("#N=",dfN[2,][[2]]), x = 12.5, y = -5)


p2 = ggplot(df, aes(x = resBs.log2FoldChange, y = -log10(resBs.padj), color = sigDir, fill = sigDir)) +
    geom_point_rast(shape = 21, stroke = 0.5, size = 1.5) +
    scale_fill_nice_udn_b() +
    scale_color_nice_udn_d() +
    geom_vline(xintercept = 0, color = "black", alpha = .5) +
    theme_nice() +
    guides(color = guide_legend(override.aes = list(size = 4))) + 
    theme(legend.key.size = unit(1, 'cm'), legend.position = "top") +
    labs(
        x = "Fold-change (log2)", 
        y = "Adjusted P value (-log10)", 
        color = "Regulation",
        fill = "Regulation") +
    # geom_text_repel(data = dfLabMa_up, aes(x = res.log2FoldChange, y = -log10(res.padj), color = sigDir, fill = sigDir, label = paste0(geneName,"-", sapply(strsplit(BsID,"\\."), `[`, 2))),
    #                 min.segment.length = 0, nudge_x = 2, size = 2.5) +
    geom_text_repel(data = dfLabMa_down, aes(x = resBs.log2FoldChange, y = -log10(resBs.padj), color = sigDir, fill = sigDir, label = paste0(geneName,"-", sapply(strsplit(BS_ID,"\\."), `[`, 2))),
                    min.segment.length = 0, nudge_x = -4, size = 2.5) +
    xlim(-6,6) +
    annotate(geom = "text", label = paste0("#N=",dfN[1,][[2]]), y = 50, x = 2.5) +
    annotate(geom = "text", label = paste0("#N=",dfN[2,][[2]]), y = 50, x = -2.5)

(p1 + p2) + plot_layout(guides = "collect") & theme(legend.position = "top")
MA and Volcano plots. With adjusted P value and fold-change cutoffs. Custom annotaitons - significant only.

Figure 10: MA and Volcano plots
With adjusted P value and fold-change cutoffs. Custom annotaitons - significant only.

# MA + Volcano plots
df = diff_bs %>%
    mutate(sig = factor(ifelse(resBs.padj < 0.05 & abs(resBs.log2FoldChange) > log2(2), TRUE, FALSE))) %>%
    mutate(dir = factor(ifelse(resBs.log2FoldChange > 0, "Up", "Down"), level = c("Up", "Down"))) %>%
    mutate(sigDir = ifelse(sig == TRUE & dir == "Up", "Up", ifelse(sig == TRUE & dir == "Down", "Down", "Not"))) %>%
    mutate(sigDir = factor(sigDir, levels = c("Not", "Up", "Down"))) %>%
    arrange(sigDir)
p1 = ggplot(df, aes(x = log2(resBs.baseMean), y = resBs.log2FoldChange, color = sigDir, fill = sigDir)) +
    geom_point_rast(shape = 21, stroke = 0.5, size = 1.5) +
    scale_fill_nice_udn_b() +
    scale_color_nice_udn_d() +
    geom_hline(yintercept = 0, color = "black", alpha = .5) +
    theme_nice() +
    guides(color = guide_legend(override.aes = list(size = 4))) + 
    theme(legend.key.size = unit(1, 'cm'), legend.position = "top") +
    labs(
        x = "Base mean", 
        y = "Fold-change (log2)", 
        color = "Regulation",
        fill = "Regulation") +
    facet_wrap(~region)

p2 = ggplot(df, aes(x = resBs.log2FoldChange, y = -log10(resBs.padj), color = sigDir, fill = sigDir)) +
    geom_point_rast(shape = 21, stroke = 0.5, size = 1.5) +
    scale_fill_nice_udn_b() +
    scale_color_nice_udn_d() +
    geom_vline(xintercept = 0, color = "black", alpha = .5) +
    theme_nice() +
    guides(color = guide_legend(override.aes = list(size = 4))) + 
    theme(legend.key.size = unit(1, 'cm'), legend.position = "top") +
    labs(
        x = "Fold-change (log2)", 
        y = "Adjusted P value (-log10)", 
        color = "Regulation",
        fill = "Regulation") +
    facet_wrap(~region)

p1 + plot_layout(guides = "collect") & theme(legend.position = "top")
MA and Volcano plots. With adjusted P value and fold-change cutoffs. Split by region.

Figure 11: MA and Volcano plots
With adjusted P value and fold-change cutoffs. Split by region.

p2 + plot_layout(guides = "collect") & theme(legend.position = "top")
MA and Volcano plots. With adjusted P value and fold-change cutoffs. Split by region.

Figure 12: MA and Volcano plots
With adjusted P value and fold-change cutoffs. Split by region.

5.1.3 Export results

# export results 
saveRDS(diff_bs, file = paste0(out, "BsDifferentialResult.rds"))

6 Session Info

sessionInfo()
## R version 4.2.2 (2022-10-31)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur ... 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] grid      stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] IHW_1.26.0                  DESeq2_1.38.3              
##  [3] SummarizedExperiment_1.28.0 MatrixGenerics_1.10.0      
##  [5] matrixStats_0.63.0          ggrastr_1.0.1              
##  [7] ggsci_2.9                   ggpointdensity_0.1.0       
##  [9] tidyr_1.3.0                 tibble_3.1.8               
## [11] patchwork_1.1.2             ggtext_0.1.2               
## [13] forcats_0.5.2               ComplexHeatmap_2.14.0      
## [15] BindingSiteFinder_1.4.0     viridis_0.6.2              
## [17] viridisLite_0.4.1           gridExtra_2.3              
## [19] ggrepel_0.9.2               kableExtra_1.3.4           
## [21] GenomicFeatures_1.50.4      UpSetR_1.4.0               
## [23] reshape2_1.4.4              dplyr_1.0.10               
## [25] AnnotationDbi_1.60.0        Biobase_2.58.0             
## [27] ggplot2_3.4.0               rtracklayer_1.58.0         
## [29] GenomicRanges_1.50.2        GenomeInfoDb_1.34.7        
## [31] IRanges_2.32.0              S4Vectors_0.36.1           
## [33] BiocGenerics_0.44.0         knitr_1.42                 
## [35] BiocStyle_2.26.0           
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.2.2               tidyselect_1.2.0         RSQLite_2.2.20          
##   [4] htmlwidgets_1.6.1        BiocParallel_1.32.5      munsell_0.5.0           
##   [7] ragg_1.2.5               codetools_0.2-18         interp_1.1-3            
##  [10] withr_2.5.0              colorspace_2.1-0         filelock_1.0.2          
##  [13] highr_0.10               rstudioapi_0.14          ggsignif_0.6.4          
##  [16] labeling_0.4.2           slam_0.1-50              lpsymphony_1.26.3       
##  [19] GenomeInfoDbData_1.2.9   mixsqp_0.3-48            polyclip_1.10-4         
##  [22] bit64_4.0.5              farver_2.1.1             vctrs_0.5.2             
##  [25] generics_0.1.3           xfun_0.36                biovizBase_1.46.0       
##  [28] BiocFileCache_2.6.0      R6_2.5.1                 doParallel_1.0.17       
##  [31] ggbeeswarm_0.7.1         clue_0.3-63              invgamma_1.1            
##  [34] locfit_1.5-9.7           AnnotationFilter_1.22.0  bitops_1.0-7            
##  [37] cachem_1.0.6             DelayedArray_0.24.0      assertthat_0.2.1        
##  [40] BiocIO_1.8.0             scales_1.2.1             nnet_7.3-18             
##  [43] beeswarm_0.4.0           gtable_0.3.1             Cairo_1.6-0             
##  [46] ensembldb_2.22.0         rlang_1.0.6              systemfonts_1.0.4       
##  [49] GlobalOptions_0.1.2      splines_4.2.2            rstatix_0.7.1           
##  [52] lazyeval_0.2.2           dichromat_2.0-0.1        broom_1.0.3             
##  [55] checkmate_2.1.0          abind_1.4-5              BiocManager_1.30.19     
##  [58] yaml_2.3.7               backports_1.4.1          Hmisc_4.7-2             
##  [61] gridtext_0.1.5           tools_4.2.2              bookdown_0.32           
##  [64] ellipsis_0.3.2           jquerylib_0.1.4          RColorBrewer_1.1-3      
##  [67] Rcpp_1.0.10              plyr_1.8.8               base64enc_0.1-3         
##  [70] progress_1.2.2           zlibbioc_1.44.0          purrr_1.0.1             
##  [73] RCurl_1.98-1.9           prettyunits_1.1.1        ggpubr_0.5.0            
##  [76] rpart_4.1.19             deldir_1.0-6             GetoptLong_1.0.5        
##  [79] ashr_2.2-54              cluster_2.1.4            magrittr_2.0.3          
##  [82] magick_2.7.3             data.table_1.14.6        circlize_0.4.15         
##  [85] truncnorm_1.0-9          SQUAREM_2021.1           ProtGenerics_1.30.0     
##  [88] hms_1.1.2                evaluate_0.20            xtable_1.8-4            
##  [91] XML_3.99-0.13            jpeg_0.1-10              shape_1.4.6             
##  [94] compiler_4.2.2           biomaRt_2.54.0           crayon_1.5.2            
##  [97] htmltools_0.5.4          Formula_1.2-4            geneplotter_1.76.0      
## [100] DBI_1.1.3                tweenr_2.0.2             dbplyr_2.3.0            
## [103] MASS_7.3-58.2            rappdirs_0.3.3           car_3.1-1               
## [106] Matrix_1.5-3             cli_3.6.0                parallel_4.2.2          
## [109] Gviz_1.42.0              pkgconfig_2.0.3          GenomicAlignments_1.34.0
## [112] foreign_0.8-84           xml2_1.3.3               foreach_1.5.2           
## [115] svglite_2.1.1            annotate_1.76.0          vipor_0.4.5             
## [118] bslib_0.4.2              webshot_0.5.4            XVector_0.38.0          
## [121] rvest_1.0.3              stringr_1.5.0            VariantAnnotation_1.44.0
## [124] digest_0.6.31            Biostrings_2.66.0        rmarkdown_2.20          
## [127] htmlTable_2.4.1          restfulr_0.0.15          curl_5.0.0              
## [130] Rsamtools_2.14.0         rjson_0.2.21             lifecycle_1.0.3         
## [133] jsonlite_1.8.4           carData_3.0-5            BSgenome_1.66.2         
## [136] fansi_1.0.4              pillar_1.8.1             lattice_0.20-45         
## [139] KEGGREST_1.38.0          fastmap_1.1.0            httr_1.4.4              
## [142] survival_3.5-0           glue_1.6.2               fdrtool_1.2.17          
## [145] png_0.1-8                iterators_1.0.14         bit_4.0.5               
## [148] ggforce_0.4.1            stringi_1.7.12           sass_0.4.5              
## [151] blob_1.2.3               textshaping_0.3.6        latticeExtra_0.6-30     
## [154] memoise_2.0.1            irlba_2.3.5.1